REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and  maintaining  the 
data  needed,  and  completing  and  reviewing  this  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggestions  for  reducing 
this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA  22202- 
4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection  of  information  if  it  does  not  display  a  currently 
valid  OMB  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 

1.  REPORT  DATE  (30-12-2011) 

2.  REPORT  TYPE  Filial 

3.  DATES  COVERED  (From  -  To) 

01  May  10-30  Apr.  11 

4.  TITLE  AND  SUBTITLE 

5a.  CONTRACT  NUMBER 

Multi-scale  modeling  of  novel  Hall  thrusters: 
understanding  physics  of  CHT  and  DCF  Thrusters 

5b.  GRANT  NUMBER 

FA9550-10-1-0188 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

Michael  Keidar 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

George  Washington  University 

8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 

9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

AFOSR,  Dr.  Mitat  Birkan 

Air  Force  Office  of  Scientific  Research 

10.  SPONSOR/MONITOR’S  ACRONYM(S) 

875  North  Randolph  St. 

Arlington,  VA  22203 

11.  SPONSOR/MONITOR’S  REPORT 

NUMBER(S) 

AFRL-OSR-VA-TR-201 2-0779 

12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Distribution  A 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

In  order  to  simulate  the  magnetized  plasma  thruster  devices  (such  as  CHT,  DCF  etc)  we  have  developed  a 
multi-scale  numerical  code  that  is  based  on  coupling  of  a  PIC/MCC  analysis  of  neutrals  and  ions  in  a 
general  2D  domain  (or  3D)  and  a  ID  kinetic  full  PIC  treatment  of  electrons  along  magnetic  fields.  The 
implementation  of  such  a  two-way  coupling  allows  calculation  of  the  electron  transport  in  the  real 
physical  domain  while  significantly  reducing  the  computational  time  associated  with  2D  full  kinetic 
simulations.  The  objective  of  this  study  was  to  identify  the  individual  contribution  to  transport  from  factors 
such  as  collisions,  surface  roughness,  secondary  electron  emission,  and  plasma  oscillations. 

Multi-scale  model  of  the  magnetized  plasma  thruster  discharge  in  which  a  kinetic  treatment  was  used  for  the  electron  component  while  a 

2D  (or  3D)  macroscopic  model  will  be  employed  for  ion  and  neutral  component  analysis.  Initial  coupling  of  the 
microscopic  and  macroscopic  model  was  performed  via  axial  electric  field,  electron  fluxes  to  the  wall 
and  electron  cross-field  transport.  This  model  demonstrated  improved  prediction  of  electron  mobility. 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF: 


17.  LIMITATION 
OF  ABSTRACT 


a.  REPORT 


c.  THIS  PAGE 


18.  NUMBER 
OF  PAGES 


19a.  NAME  OF  RESPONSIBLE  PERSON 


19b.  TELEPHONE  NUMBER  (include  area 
code) 


b.  ABSTRACT 


Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  Z39.18 


Multi-scale  modeling  of  novel  Hall  thrusters:  understanding  physics  of  CHT  and 


DCF  thrusters 


Michael  Keidar 


Department  of  Mechanical  and  Aerospace  Engineering 
The  George  Washington  University,  Washington  DC 


Final  Report 


AFOSR  Grant  Number:  FA9550- 10-1-0188 


Abstract 


In  order  to  simulate  the  magnetized  plasma  thruster  devices  (such  as  CHT,  DCF  etc)  we  have  developed  a 
multi-scale  numerical  code  that  is  based  on  coupling  of  a  PIC/MCC  analysis  of  neutrals  and  ions  in  a 
general  2D  domain  (or  3D)  and  a  ID  kinetic  full  PIC  treatment  of  electrons  along  magnetic  fields.  The 
implementation  of  such  a  two-way  coupling  allows  calculation  of  the  electron  transport  in  the  real 
physical  domain  while  significantly  reducing  the  computational  time  associated  with  2D  full  kinetic 

simulations.  The  objective  of  this  study  was  to  identify  the  individual  contribution  to  transport  from 
factors  such  as  collisions,  surface  roughness,  secondary  electron  emission,  and  plasma  oscillations.  Multi¬ 
scale  model  of  the  magnetized  plasma  thruster  discharge  in  which  a  kinetic  treatment  was  used  for  the 
electron  component  while  a  2D  (or  3D)  macroscopic  model  will  be  employed  for  ion  and  neutral 
component  analysis.  Initial  coupling  of  the  microscopic  and  macroscopic  model  was  performed  via  axial 
electric  field,  electron  fluxes  to  the  wall  and  electron  cross-field  transport.  This  model  demonstrated 
improved  prediction  of  electron  mobility. 


Section  1 ,  Multiscale  Modeling  of  Hall  Thrusters 


Introduction 

Despite  Hall  thrusters  having  over  40  years  of  flight  heritage  (the  first  variant,  SPT-50,  was  flown 
aboard  the  Soviet  Meteor  spacecraft  in  1971),  the  community  still  lacks  a  tool  capable  of  predictively 
modeling  these  devices.  The  closest  to  an  industry  standard,  COTS-level  simulation  program  is 
HPHall.  This  2D  axisymmetric  code  was  originally  developed  at  MIT  by  Fife,1  and  was  subsequently 
improved  by  others.2’ 3  Other  researchers  have  developed  their  own  codes  based  on  the  basic  premise 
of  HPHall.4-7  These  codes  share  a  common  approach.  Ions  are  treated  as  kinetic  particles  while 
electrons  are  modeled  as  a  quasi- ID  fluid.  Quasineutrality  is  assumed,  and  the  kinetically 
determined  plasma  density  is  used  to  update  the  electron  temperature  in  direction  perpendicular  to 
the  magnetic  field  lines.  Potential  is  recovered  by  considering  current  conservation  and  assuming 
electron  thermalization  along  field  lines.  This  model  is  generally  valid  -  in  Hall  thrusters  electrons 
are  magnetized,  and  thus  their  motion  can  be  decoupled  into  two  distinct  modes:  the  unhindered 
motion  along  the  field  line,  and  the  mobility-driven  diffusion  across  them.  The  mobility  term  comes 
into  play  in  the  electron  momentum  equation  used  to  determine  the  cross-field  velocity.  This  is  one  of 
the  terms  in  the  energy  equation.  HPHall  relies  on  the  classical  1/B2  collision  driven  transport 
combined  with  an  anomalous  1  /B  Bohm  term.  Mathematically,  the  mobility  model  is  ge,n  =  pe  /(1+ 
P2)  +  Kb  /(16B),  where  \iQis  the  mobility  \iQ=  e/venmQ  for  unmagnetized  electrons,  Kb  is  the  user- 
determined  Bohm  parameter  and  |3  =  <dc  /v  is  the  Hall  parameter.  This  formulation  is  not  without 
problems.  Transport  is  inherently  a  kinetic  phenomenon.  It  is  influenced  to  a  large  degree  by  wall 
collisions,  secondary  electron  emissions,  and  temporal  field  variations.  These  are  kinetic  effects  that  the 
simple  analytical  model  fails  to  capture.  In  addition,  as  pointed  out  previously  by  Sydorenko,8 
electron  temperature  in  Hall  thrusters  is  anisotropic,  further  complicating  the  model. 

Hence,  several  other  researchers  have  taken  path  completely  opposite  of  HPHall,  and  have  developed 
fully  kinetic  Hall  thruster  codes.  Example  of  one  such  code  is  the  work  of  Hirakawa.9  Many  other 
researchers  have  also  contributed.10-12  These  fully  kinetic  codes  have  their  own  drawbacks.  First,  they 
require  extremely  small  simulation  time  steps.  Resolving  the  cyclotron  rotation  requires  timesteps 
some  six  orders  of  magnitude  smaller  than  required  to  push  ions  in  HPHall.  Second,  since  potential  in 
these  codes  is  obtained  from  solution  of  the  Poisson’s  equation,  the  computational  cells  must  be 
small  enough  to  resolve  the  Debye  length.  These  two  requirements  make  fully  kinetic  approach  a 
computationally  daunting  task  for  all  but  the  smallest  of  thrusters.  Even  if  tricks  are  played  with 
relative  electron-to-ion  mass  and  the  free  space  permittivity,  these  codes  still  require 
supercomputers  and  many  days  or  weeks  of  computational  effort. 

Fully  kinetic  codes  are  useful  tools  for  understanding  basic  physical  processes  occurring  in  Hall 
thrusters.  They  are  not  particularly  useful,  however,  to  a  designer  working  on  optimizing  one  of 
these  devices.  Designer  at  some  aerospace  corporation  will  likely  not  have  access  to  supercomputing 
facilities,  nor  have  the  time  to  perform  a  trade  study  if  each  cases  requires  substantial  computational 
effort.  For  this  reason,  we  started  working  on  a  new  multiscale  approach  to  modeling  Hall  thrusters. 
The  objective  of  our  work  is  to  develop  a  tool  capable  of  self-consistently  determining  electron 
mobility  and  the  thruster  plasma  properties  of  interest,  but  do  it  in  such  a  way  that  a  solution  can  be 
obtained  on  a  standard  desktop  workstation  in  less  than  a  day.  To  accomplish  this  goal,  we  divide 
the  problem  into  the  following  three  spatial  scales: 

•  Magnetic  Field  Line:  On  the  spatial  scale  of  a  magnetic  field  lines,  dynamics  is  driven  by 
the  cyclotron  motion  of  electrons.  Electrons  are  magnetized,  and  individual  field  lines  can  be 


considered  independent  of  each  other.  Heavy  particle  and  properties  normal  to  the  field  lines  are 
assumed  frozen.  This  approach  allows  us  to  rapidly  simulate  electrons  and  recover  mobility  self- 
consistently.  Leveraging  modern  multi-  core  architectures  via  multithreading  allows  us  to  study 
multiple  field  lines  simultaneously. 

•  Thruster  Channel:  On  the  spatial  scale  of  the  thruster,  plasma  is  assumed  to  be 
quasineutral,  and  electron  density  can  be  obtained  from  kinetic  ions.  Electron  temperature  and 
plasma  potential  is  obtained  by  solving  the  quasi- ID  equations.  This  approach  is  identical  to 
HPHall,  except  that  our  method  relies  on  the  mobility  determined  by  the  kinetic  analysis.  Ions 
exiting  the  thruster  at  steady  state  are  sampled  to  obtain  a  discretized  source  term  for  plume 
modeling. 

•  Plume  Environment:  Outside  the  thruster  exit,  the  magnetic  field  plays  a  negligible  role. 
The  plume  is  quasineutral  except  in  low  density  sheath  regions  around  the  spacecraft.  Of  interest 
here  is  the  formation  of  charge  exchange  ions  and  their  impact  on  spacecraft  components.  Electron 
density  is  obtained  from  Boltzmann  relationship,  and  potential  can  be  solved  by  direct  inversion 
or  by  solving  Poisson’s  equation. 


lines  studied  in  Step  1 
are  shown  in  red. 

A/ 

Figure  1.  Schematic  of  our  multiscale  approach.  Codes  1.  and  2.  are  designed  to  iterate  until  convergence.  We  then 
sample  ions  crossing  the  thruster  exit  plane  to  obtain  the  source  model  for  plume  analysis. 

This  approach  is  shown  schematically  in  Figure  1.  This  formulation  naturally  lends  itself  to  three 
codes,  each  concentrating  on  the  physics  of  the  respective  spatial  scale.  The  first  code,  Lynx, 
simulates  the  cyclotron  motion  of  electrons  about  series  of  magnetic  field  lines.  The  second  code,  a 
general  2D  axisymmetric  /  Cartesian  plasma  simulation  tool  Starfish,  is  still  in  development,  and 
hence  in  this  paper  we  utilized  HPHall.  The  final  piece  is  the  3D  plume  code  Draco.  We  begin  the 
paper  by  describing  the  thruster  to  which  this  approach  was  applied.  We  next  skip  to  the  second  step, 
and  use  HPHall  to  obtain  the  initial  plasma  parameters  which  will  serve  as  inputs  to  Lynx.  We  then 
perform  a  kinetic  analysis  to  compute  new  values  mobility.  These  are  then  used  to  update  the  HPHall 
solution  and  extract  the  plume  source  model.  Repetitive  iteration  between  steps  1  and  2  is  left  as  a 
future  item  that  will  be  attempted  once  Starfish  development  is  complete.  The  paper  concludes  with 
a  summary  and  a  discussion  of  future  work. 


Thruster  Details 


We  deploy  our  model  to  the  2.6cm  Princeton  Cylindrical  Hall  thruster  (CHT).  This  thruster  is 
described  in  greater  detail  in  Ref.  13.  Here  we  summarize  just  the  parameters  important  to  our 
analysis.  The  most  important  characteristic  of  this  device  is  its  non-standard  geometry.  While  typical 
Hall  thrusters  consist  of  an  annular  channel,  the  CHT  contains  an  annular  upstream  zone  and  a 
cylindrical  acceleration  zone.  The  lack  of  the  inner  wall  in  the  acceleration  region  is  expected  to  lead 
to  an  increased  thruster  lifetime  and  improved  performance  due  to  reduced  losses  of  ions  to  the 
walls.  From  the  academic  standpoint,  this  configuration  also  introduces  interesting  new  physics.  The 
magnetic  field  lines  converge  near  the  innerpole,  resulting  in  a  region  of  increased  magnetic  pressure. 
Electrons  are  then  expected  to  be  preferentially  scattered  to  the  outer  wall,  possibly  resulting  in  a 
non-symmetric  sheath.  This  finding  was  touched  upon  in  our  previous  work.14 

The  walls  of  this  SPT-type  thruster  are  made  of  a  dielectric  material.  The  dielectric  walls 
distinguish  this  thruster  design  from  another  Hall  thruster  variant,  TAL  (thruster  with  anode 
layer),  in  which  the  walls  are  conductive.  One  theory  of  SPT  operation  suggests  that  electron  impacts 
of  the  dielectric  material  result  in  emission  of  secondary  electrons  from  the  material  matrix.  Since 
these  electrons  are  not  initially  magnetized,  they  are  free  to  migrate  towards  the  anode.  This  theory  is 
called  near  wall  conductivity  (NWC)  and  it’s  a  feature  that  is  captured  by  our  kinetic  code.  These 
electrons  are  also  significantly  colder  than  the  impacting  particles,  resulting  in  the  cooling  of  the 
primary  population. 

Experimental  measurements  of  this  thruster  have  been  presented  in  Refs.  15  and  16.  The  CHT  can 
be  operated  in  two  modes  based  on  the  current  applied  to  the  magnetic  coils.  Our  work  correlates  to 
the  ’’direct”  mode.  In  this  configuration,  magnetic  field  lines  cross  the  channel  between  the  outer  wall 
and  the  innnerpole  without  forming  a  cusp  near  the  outer  wall.  Electron  temperature  in  the  thruster 
was  found  to  peak  at  25  eV  just  outside  the  exit  plane.  Potential  decays  slowly  through  majority  of 
the  channel.  Only  50V  of  potential  drop  were  measured  in  the  first  1.5cm  from  the  anode.  Most  of 
the  potential  drop  occurs  near  the  exit  plane,  in  the  acceleration  zone.  Additional  50V  potential 
drop  was  measured  to  occur  outside  the  thruster.  Anode  current  was  approximately  0.3A,  and 
plasma  density  was  ~6X1017m3 . 

Initial  Results  for  Thruster  Discharge 

We  used  the  HPHall  simulation  code  to  model  the  discharge  channel.  Schematic  of  our  setup  is  shown 
in  Figure  2.  HPHall  uses  a  structured  but  non-uniform  mesh,  boundary  of  which  is  shown  in  the 
figure.  Such  a  mesh  simplifies  capturing  the  geometry  of  non-rectangular  devices,  as  well  as  the 
downstream  near-plume  region.  Internally,  HPHall  establishes  another  virtual  mesh  that  is  used  to 
solve  the  electron  equations.  This  mesh  is  illustrated  in  Figure  2.  The  vertical  lines  correspond  to 
approximately  equidistantly  spaced  magnetic  field  X  lines.  It  should  be  noted  that  this  graphics  does 
not  completely  correlate  to  the  implementation  in  HPHall.  The  data  structure  used  internally  by  the 
code  uses  a  variable  number  of  radial  segments,  which  leads  to  a  more  complicated  visualization 
problem.  The  mesh  shown  here,  and  in  our  subsequent  analysis,  used  a  uniform  number  of  radial 
segments  corresponding  to  the  average  number  used  by  HPHall. 

Of  importance  are  the  left  and  the  right  boundaries.  These  correspond  to  the  anode  and  the  cathode, 
respectively.  Constant  potential  is  applied  upstream  of  the  anode.  The  anode  potential  is  computed 
by  the  code  self  consistently.  The  right  boundary  is  the  cathode  line  on  which  potential,  density,  and 
temperature  are  specified.  Simple  linear  interpolation  for  temperature  and  potential  is  used 
downstream  of  the  cathode.  As  can  be  seen  from  the  figure,  this  downstream  region  can  encompass 
a  substantial  fraction  of  the  simulation  domain.  Important  part  of  setting  up  an  HPHall  simulation  is 


determining  where  to  place  the  anode  and  the  cathode  lines.  Typically,  trade  study  is  performed  and 
the  solution  in  the  best  agreement  with  experiments  is  selected. 


Figure  2.  Simulation  model  of  the  cylindrical  Hall  thruster.  Slice  shows  the  HPHall  computational  domain.  The 
mesh  corresponds  to  the  region  in  which  the  electron  conservation  equations  are  solved. 


HPHall  injects  kinetic  electrons  at  the  anode.  Ionization  model  is  then  used  to  create  ion  particles.  Ion 
positions  are  integrated  using  the  particle  in  cell  method.  The  code  assumes  charge  neutrality,  and  hence 
71  e  ~  ni  +  2nl2  ,  summation  of  singly  and  doubly  charged  ion  populations.  In  order  for  quasineutrality 

to  hold,  electron  temperatures  must  adjust  accordingly  to  accommodate  diffusion.  HPHall  solves  the 
energy  equation 


d  /3 

di 


G  kT° )  +  V  +  Qe)  +  v  '  ( nekTeue )  =  Sh-Si 


(1) 


on  the  previously  described  lambda  mesh.  Constant  electron  temperature  at  each  field  line  is  assumed. 
This  then  allows  each  field  line  to  be  treated  as  a  single  volume  element.  Relevant  properties,  such  as 
mobilities  or  electron  densities,  are  integrated  along  the  line  following  the  standard  finite  volume 
formulation.  In  other  words,  the  computation  is  performed  in  a  quasi- ID  dimension.  The  conservation 
equations  are  solved  only  in  the  direction  normal  to  the  field  lines,  but  2D  radial  contribution  is  used  to 
compute  the  coefficients  at  each  point.  Electron  velocity,  given  by  the  momentum  balance, 


Ue  =Fe 


(2) 


is  incorporated  into  the  temperature  solver.  Here  pe  is  the  mobility  term.  Once  temperature  has  been 

determined,  the  thermalized  potential  cp  can  be  computed  at  each  lambda  line  from  current  conservation. 
The  radial  variation  in  potential  is  then  recovered  from  the  thermalized  potential  relationship, 
(p-  =  <p-  kTe/eln(ne). 


In  our  previous  work14  we  performed  a  limited  parametric  study  where  we  investigate  several  different 
cathode  line  positions  and  Bohm  parameters  with  the  goal  of  matching  the  experimental  measurements 
for  CHT.  We  were  able  to  improve  our  correlation  for  this  work.  The  relevant  input  settings  are 
summarized  in  Table  1.  Plasma  properties  along  the  thruster  centerline  are  plotted  in  Figure  6.  This  figure 
shows  both  the  results  obtained  using  the  analytical  mobility  model,  and  the  results  obtained  subsequently 
using  our  kinetic  method.  The  initial  results  with  the  analytical  mobility  are  plotted  with  dashed  lines.  It 
should  be  noted  that  there  are  still  discrepancies.  For  one,  the  HPHall-predicted  potential  distribution 
fails  to  capture  the  gradual  rise  to  the  anode  potential.  Instead,  potential  is  actually  seen  to  decrease 
upstream  of  the  start  of  the  acceleration  zone.  This  indicates  that  fraction  of  ions  born  in  this  regions  will 
have  the  tendency  to  flow  towards  the  anode.  One  of  the  goals  for  this  work  was  to  determine  if  any 
improvement  can  be  achieved  by  modifying  the  mobility  distribution  using  the  kinetic  code. 

Kinetic  Code  Inputs 

The  kinetic  code  requires  as  inputs  information  related  to  the  global  state  of  the  discharge.  We  developed 
a  simple  code  to  contour  HPHall  results  (2d_ave_tp.dat)  along  the  lambda  lines  for  computation. 
Contouring  starts  by  searching  for  the  corresponding  value  along  the  bottom  edge  of  the  computational 
domain.  Edge  cuts  are  then  determined  by  linear  interpolation  of  node  values,  and  the  cuts  are  connected 
to  form  a  spline.  Properties  of  interest  are  then  interpolated  onto  the  spline  control  points.  Figure  3  shows 
the  parameters.  Along  with  the  magnetic  field  profile  (not  shown  here),  these  five  parameters  serve  as 
inputs  to  the  kinetic  code.  The  kinetic  code  is  described  in  the  next  section. 

Table  1.  Summary  of  critical  inputs  for  the  HPHall  simulation 


Parameter 

Value 

Mass  flow  rate 

4  x  10-7  (kg/s) 

Discharge  voltage 

275  V 

Anode  line  position 

(0.0520,  0.010)  m 

Cathode  line  position 

(0.0670,  0.010)  m 

Cathode  temperature 

26  eV 

Cathode  potential 

210.6  V 

Cathode  density 

xl017  m-3 

Cathode  emitter  potential 

-20  V 

Bohm  coefficient 

1 

(a)  Plasma  Density 


(b)  Neutral  Density 


(c)  Normal  Component  of  Electric  Field 


(d)  Electron  Temperature 


Figure  3.  Plasma  parameters  serving  as  inputs  to  the  kinetic  code.  Kinetic  simulation  is  performed  for  each  field 
line  (vertical  grid  lines).  Additional  input,  which  is  not  shown  here,  is  the  magnetic  field  strength. 


Mobility  Calculation  with  Lynx 

The  inputs  in  Figure  3  are  next  processed  with  a  kinetic  code  called  Lynx.  Lynx  is  a  kinetic  code  that  self- 
consistently  calculates  the  radial  variation  in  electron  mobility  along  a  magnetic  field  line.14,  17  Lynx 
simulates  only  electrons.  Ions  and  neutrals  are  assumed  to  remain  frozen  during  the  time  necessary  to 
compute  electron  trajectories.  The  heavy  particles  thus  form  a  fixed  background  with  which  the  electrons 
interact  during  their  cyclotron  motion  about  the  field  line.  By  repeating  the  calculation  for  all  the  lambda 
lines  making  up  the  HPHall  electron  mesh,  we  can  obtain  a  two-dimensional  variation  in  mobility.  This 
kinetically  determined  mobility  can  then  be  used  in  HPHall  in  lieu  of  the  analytical  model. 

Lynx  is  implemented  in  Java.  Although  historically  Java  performance  was  not  competitive  with  languages 
such  as  C/C++,  this  is  no  longer  the  case.  Modem  Java  compilers  generate  codes  that  perform  at  speeds 
comparable  and  in  some  cases  even  exceeding  native  C++  implementation.  1 8  In  addition,  Java  offers  a 
large  standard  library  with  support  for  data  management,  as  well  as  GUI,  networking,  and  graphics 
rendering.  In  addition,  Java  natively  supports  multithreading  which  allows  us  to  take  advantage  of  modern 
multi-core  system  architectures,  and  mn  Lynx  concurrently  for  multiple  field  lines.19  Similar  functionality 
can  be  achieved  in  C++  through  the  Boost  libraries.  However,  that  step  requires  downloading, 
configuring,  and  compiling  the  massive  libraries.  In  Java,  such  support  is  provided  natively. 

A.  Execution  and  Topology 

Lynx  simulation  commences  with  the  code  importing  the  2D  mesh  shown  previously  in  Figure  3.  An 
additional  settings  file  is  also  processed.  This  settings  file  specifies  general  simulation  parameters  such  as 
wall  temperature,  number  of  electrons  per  field  line,  and  the  number  of  time  steps.  It  also  controls  which 
diffusion-inducing  processes,  such  as  collisions  or  wall  effects  should  be  included.  Lynx  next  instantiates 
a  new  simulation  object  for  each  magnetic  field  line.  A  simple  scheduler  was  implemented  to  launch  the 
threads  and  monitor  them  for  completion.  The  number  of  threads  running  concurrently  is  limited  by  the 
number  of  available  CPU  cores. 


The  computational  domain  for  each  Lynx  simulation  thread  is  a  one  dimensional  domain  corresponding  to 
a  particular  lambda  line.  In  Figure  3,  each  vertical  grid  line  corresponds  to  an  individual  simulation.  The 
number  of  nodes  along  the  line  differed  from  the  coarse  mesh  shown  in  this  figure  and  was  determined 
automatically  such  that  As  =  0.5?JT  .  On  average,  100  nodes  were  used  per  magnetic  field  line.  Note,  the 
coarse  input  mesh  was  used  to  store  transport  properties.  The  coarser  mesh  was  selected  in  order  to  reduce 
statistical  noise  errors.  Several  fixed  properties  are  set  at  each  grid  node  during  the  initialization.  These 


include  ni  and  nO  ,  ion  and  neutral  densities,  E  ,  perpendicular  component  of  electric  field,  B,  magnetic 
field  strength,  and  SB/5s,  magnetic  field  gradient. 

B.  Particle  Loading 

The  simulation  then  continues  by  loading  the  electron  particles.  Electrons  are  loaded  at  grid  nodes,  with 
the  number  of  electrons  obtained  by  multiplying  the  ion  density  by  cell  volume  and  scaling  by  the  particle 
weight.  As  pointed  out  by  Fox,  11  careful  sampling  of  electrons  is  necessary  in  order  to  resolve  the  high 
energy  tail  of  the  velocity  distribution  function.  We  did  not  take  this  approach  at  this  time,  but  plan  to 
investigate  the  effect  of  variable  particle  weight  on  transport  and  near  wall  conductivity  in  the  near  future. 
The  initial  electron  temperature  at  each  field  line  is  one  of  the  inputs  from  HPHall.  Isotropic  temperature 
distribution  was  used  for  particle  loading.  Particles  are  loaded  such  that  the  axial  position  of  the  guiding 
center  is  z  =  0.  Initial  number  of  particles  was  50,000. 

The  initial  plasma  sheath  profile  needs  to  be  determined  next.  Since  ions  are  frozen,  the  establishment  of 
equilibrium  sheath  is  marked  by  zero  electron  current  to  the  walls.  During  this  initial  period,  collisions 
were  not  performed,  and  secondary  electrons  were  not  emitted.  Instead,  particles  impacting  the  walls  were 
simply  removed  from  the  simulation  domain.  These  processes  were  ignored  since  we  are  interested  in 
determining  the  electron  density  distribution  arising  from  the  thermal  and  potential  balance.  In  our 
simulations,  we  found  that  steady  state  was  achieved  after  approximately  1500  time  steps.  However,  for 
an  added  margin  of  safety,  we  required  to  code  to  complete  3000  steps  before  the  simulation  continued. 
The  additional  time  steps  have  no  impact  on  the  results  since  once  the  steady  state  is  achieved,  electrons 
become  trapped  in  the  potential  well  between  the  walls  and  are  unable  to  reach  the  walls. 

Simulation  then  moves  to  the  normal  operating  mode,  with  collisions  and  SEE  turned  on.  Particles 
translating  away  from  the  field  line,  or  impacting  the  walls  were  removed  from  the  simulation.  The  code 
reinjected  new  particles  into  the  simulation  to  keep  the  number  of  electrons  at  the  steady  state  level.  Total 
of  20,000  time  steps  were  simulated.  Averaging  of  results  began  at  time  step  5000.  The  intermediate 
period  between  steady  state  and  averaging  was  excluded  from  averaging  to  avoid  any  possible  initial 
transient  effects  from  influencing  the  results.  The  timestep  At  was  selected  such  that  electrons  completed 
single  orbit  in  75  time  steps. 

C.  Collisions  and  Near  Wall  Conductivity 

Electron  collisions  were  modeled  using  the  Monte  Carlo  method.  In  this  method,  source  particles  are 
collided  with  a  stationary  target  cloud.  The  collision  probability  is  determined  from  the  background 
density,  nO  and  collision  frequency  is  P  =  1  —  exp(—n0  a0  gAt )  Here  a0  is  the  total  collision  cross- 

section  due  to  all  processes.  For  particles  undergoing  collision,  the  collision  process  was  picked  randomly 
according  to  the  ratio  of  ax  /o0.  Post-collision  velocity  was  computed  by  first  sampling  a  random  target 

velocity  and  a  random  impact  angle.  We  then  calculated  the  post-collision  velocity  from  conservation  of 
energy. 

To  reduce  statistical  errors  and  improve  performance,  collisions  were  computed  only  once  every  4 
timesteps.  Three  types  of  collisions  were  considered:  momentum-transfer  (electron-atom),  ionization 
(electron-atom),  and  excitation  (electron-atom).  Coulomb  collisions  were  not  included  in  the  present 
work.  At  low  electron  temperatures,  polarization  collisions  dominate  the  momentum-transfer  interaction 
between  electrons  and  atoms.20  Cross-section  for  this  process  was  obtained  from  the  analytical  model20 
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where  ar  is  the  polarizability  of  the  atom.  It  is  given  by  ar  =  27.66a3  ,  where  a0  is  the  Bohr  radius.21 
Cross-section  for  electron-ion  collisions  was  given  by 
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where  b0  =  q^q2  /( 2ne0  mv2  )  is  the  distance  of  closest  approach.  The  inelastic  process,  ionization  and 

ex-  citation,  cross-sections  were  computed  using  the  polynomial  fit  of  Szabo.22  These  collisions  were 
modeled  by  reducing  the  energy  of  the  impacting  electron  by  the  ionization  or  excitation  energy,  and 
scattering  the  electron  through  a  random  angle. 


Electrons  impacting  the  dielectric  walls  were  reflected  back  to  the  domain,  absorbed,  or  generated  sec¬ 
ondary  electrons  according  to  the  model  of  Sydorenko.23  The  secondary  electron  yield  was  given  by, 


Y  = 


(5) 


The  SEE  yield  is  non-negligible  for  kTe  =  lOeV.  The  SEE  electrons  were  assumed  to  come  off  the  surface 
unmagnetized,  and  were  generated  at  the  wall  with  initial  direction  given  by  a  random  velocity  vector. 
Impacting  electron  knocking  off  a  secondary  electron  was  assumed  to  be  absorbed  by  the  wall  to  retain 
charge  neutrality,  and  was  removed  from  the  simulation.  Although  we  have  implemented  a  simple 
analytical  model  to  take  into  account  surface  charging  of  dielectric  walls,  we  did  not  utilize  it  at  present. 
Instead,  potential  was  fixed  at  both  walls  at  OV. 

Figure  4  shows  the  typical  field  line  quantities.  This  particular  plot  was  generated  for  a  magnetic  field  line 
close  to  the  anode.  We  can  see  that  in  the  bulk  region,  electrons  clearly  follow  the  ions.  This  charge 
neutrality  is  captured  by  the  potential  solution,  leading  to  zero  electric  field  in  this  bulk  region.  A  clear 
sheath  forms  near  the  walls.  Electron  density  decays  faster,  as  expected.  Electric  field  forms  near  the 
walls  to  repel  low  energy  electrons.  Wall  potential  drop  is  5.7kTe. 

D.  Transport  Calculation 

Previous  paragraphs  described  the  methods  used  to  simulate  magnetized  electrons  bounded  between  two 
walls.  In  order  to  make  the  model  useful  to  our  multiscale  approach,  we  need  to  extract  mobility  from  the 
solution.  Conceptually,  determining  mobility  is  a  trivial  task.  Mobility  could  be  obtained  as 
Hj  =  ZUz/EL  where  j  corresponds  to  a  mesh  cell,  and  the  sum  is  performed  over  all  particles  in  that 

cell.  Unfortunately,  this  approach  is  not  feasible  numerically.  To  see  why,  let’s  consider  the  typical 
values  of  drift  velocity  vd  and  the  tangential  velocity  of  electrons  orbiting  a  field  line,  vT  .  Using 
values  typical  of  the  CHT,  vd  =  pE  =  5m2  /Vs  X  20,000V/m  =  105m/s.  The  tangential  velocity  can  be 


approximated  as  vT  —  ^JkTe/m  =  106m/s.  From  this  simple  calculation,  we  can  see  that  vT  >  vd. 

Given  a  finite  number  of  particles,  and  the  fact  the  particles  translate  along  the  field  line  during  their 
orbit,  the  summation  will  always  result  in  non-zero  mobility.  Due  to  the  magnitude  of  vT  ,  a  single 
unpaired  term  will  likely  wash  out  any  actual  drift  velocity. 

This  shortcoming  became  obvious  during  our  validation  tests.  We  ran  the  simulation  for  a  case  in 
which  all  diffusion  terms  were  disabled,  yet  the  code  predicted  finite  transport.  Even  more 
interesting  was  the  fact  that  mobility  was  only  slightly  increased  when  collisions  were  included  or 
when  electric  field  was  doubled.  Hence,  we  implemented  an  alternative  method  of  computing  drift 
velocity.  Instead  of  summing  the  axial  components  of  particle  velocity,  we  consider  the  velocity  of  the 
guiding  center.  We  determine  the  guiding  center  as  rg  =  0.5(z+  +  z~),  where  z+  and  z-  are  the 
extents  of  particle’s  position  during  an  orbit.  These  quantities  are  reset  once  per  orbit.  Drift  velocity 
is  then  vd=  rg/x,  where  x  is  the  time  delta  since  previous  sampling. 


Figure  4.  Plot  of  typical  simulation  results.  Electrons  are  seen  to  closely  follow  the  prescribed  ion  density ,  except  in 
the  sheath  region,  where  they  density  decays  more  rapidly  as  expected.  Plasma  potential  and  tangential  electric 

field  profiles  are  also  shown. 

Only  particles  having  diffused  more  than  1.5rL  were  counted.  Here  rL  is  the  local  Larmor  radius 
computed  using  particle’s  tangential  velocity,  and  the  strength  of  the  magnetic  field  at  the  particle 
location.  Stray  particles  were  removed  from  the  simulation  and  were  subsequently  replaced  by  new 
particles  at  the  field  line.  The  radial  position  of  the  newly  created  particles  was  based  on  the  density 
of  the  ions. 

All  simulations  presented  in  this  paper  were  executed  on  a  Sony  VAIO  laptop  with  the  Intel  Core  i7 
2.5GHz  CPU.  This  CPU  supports  up  to  8  concurrent  threads.  Typical  simulation  times  ranged  from 
400  seconds  for  a  quick  estimate  with  10,000  particles  and  10,000  time  steps  to  2600  seconds  for  the 
mns  used  to  generate  the  results  for  this  paper.  These  runs  used  50,000  particles  and  ran  for  20,000 
time  steps.  Computed  properties,  including  mobility,  were  then  exported  as  a  2D  mesh. 

Figure  5  compares  the  mobility  computed  using  our  kinetic  approach  to  that  used  by  HPHall. 
Although  the  background  values  of  mobility  are  similar  quantitatively,  we  can  see  stark  differences 
between  the  two  versions.  The  kinetic  solution  contains  two  distinct  regions  of  high  mobility  which  are 
not  seen  in  the  analytical  model.  High  production  of  secondary  electrons  was  predicted  by  the  code 
for  the  field  lines  at  the  left  band.  NWC  may  also  explain  the  oscillatory  nature  of  transport,  which 
seems  to  be  related  to  the  high  number  of  SEE  seen  on  the  magnetic  field  lines  in  this  region.  On  the 


other  hand,  the  high  mobility  in  the  right  band  may  be  due  to  a  strong  electric  field.  The  right  field 
line  corresponds  to  the  location  where  HPHall  predicts  drop  in  potential  corresponding  the  the  start 
of  the  acceleration  zone.  We  also  see  reduced  mobility  near  the  innerpole.  This  region  is  dominated 
by  increased  magnetic  pressure  which  reduces  flux  of  electrons  to  this  region.  It  should  be  noted  that 
the  results  shown  here  are  statistically  accurate.  We  performed  multiple  simulations  with  a  varying 
number  of  particles,  and  all  simulations  produced  comparable  results. 


(a)  HPHall  (b)  Kinetic  Code 

Figure  5.  Comparison  of  the  analytical  mobility  used  by  HPHall  (left)  to  the  kinetically  determined  mobility  from 

Lynx. 

HPHall  results  with  kinetic  mobility 

The  mobility  contour  obtained  in  Lynx  was  next  loaded  into  HPHall.  We  modified  HPHall  to  load 
mobility  values  from  a  file.  We  next  ran  HPHall  for  5,000  time  steps,  and  exported  the  new  2D  results. 
Comparison  of  centerline  plasma  properties  can  be  seen  in  Figure  6.  Solid  markers  correspond  to 
experimental  data.  The  dashed  lines  are  the  results  obtained  with  the  classical  model,  and  the  solid 
lines  are  the  ones  obtained  by  loading  the  self-consistently  determined  mobility.  A  clear 
improvement  in  plasma  potential  is  seen.  Although  the  potential  deviates  somewhat  from  the  control 
point,  the  potential  is  seen  to  increase  towards  the  anode,  and  the  trough  region  is  eliminated. 
Temperature  results  are  less  conclusive,  partly  due  to  the  poor  agreement  of  the  initial  results  with 
the  data.  The  peak  of  plasma  density  is  seen  to  decrease  and  move  towards  the  anode. 
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Figure  6.  Plot  of  several  plasma  parameters  along  the  thruster  centerline.  The  markers  correspond  to  experi¬ 
mental  data  from  Ref  15  Experimental  plasma  density,  measured  in  a  single  point  in  the  channel,  was  6  x  1017 
m~3.  The  results  computed  using  the  analytical  mobility  model  are  shown  using  the  dashed  lines.  The  results 
obtained  with  the  kinetic  mobility  are  shown  with  the  solid  lines. 


A.  Particle  Sampling 

During  this  simulation  we  also  sampled  particles  crossing  a  virtual  plane  in  the  near  plume  region.  It’s 
difficult  to  describe  the  velocity  and  flux  distribution  function  of  the  Hall  thruster  plume  analytically. 
Unlike  ion  thrusters,  Hall  thruster  discharges  are  open  to  the  ambient  environment,  and  the 
acceleration  profile  is  a  function  of  the  plume  itself.  Sampling  particles  as  they  cross  some  plane  is  a 
simple  and  effective  way  to  describe  the  velocity  and  flux  space.  By  sampling  velocities  and  the 
corresponding  spatial  position  for  asufficiently  large  number  of  particles  we  can  obtain  a  discretized 
velocity  distribution  function.  A  peculiar  feature  of  Hall  thrusters  is  that  some  acceleration  occurs 
outside  the  actual  thruster.  In  the  case  of  the  CHT,  the  external  acceleration  accounts  for 
approximately  50eV  of  ion  energy.  In  order  to  capture  this  feature,  we  sample  the  particles  sufficient 
distance  from  the  exit  plane,  in  a  region  where  potential  drop  becomes  negligible. 

Plume  simulations  with  Draco 

The  discretized  velocity  distribution  function  is  next  used  to  model  the  plume  environment  produced 
by  these  thrusters.  The  plume  modeling  was  performed  using  the  final  piece  of  our  multiscale 
approach,  a  3D  ES-PIC  code  called  Draco.24  Draco  was  developed  at  Virginia  Tech  as  a  general 
plasma  simulation  tool  suitable  for  modeling  electric  propulsion  plumes  and  their  interaction  with 
spacecraft  components.  It  readily  interfaces  with  detailed  geometries  produced  in  COTS  CAD/CAE 
packages.  The  code  operates  on  a  rectilinear  mesh  containing  a  subset  of  cut  cells.  The  cut  cells  are 
generated  automatically  based  on  the  user  specified  geometry  surface  mesh. 

A.  Source  Model 

Particles  are  injected  into  the  simulation  from  source  element  groups.  Draco  supports  variety  of  source 
models,  including  one  developed  specifically  for  interfacing  with  HPHall.25  This  model  takes  as  its 
input  a  file  containing  a  large  number  of  [r,  vz ,  vr,  v0]  tuples.  The  source  randomly  selects  a  tuple  from 
the  list  and  rotates  it  through  a  random  azimuthal  angle  0.  Position  and  and  velocity  is  also  rotated 
according  to  the  normal  vector  of  the  source  elements.  This  discretized  model  offers  the  benefit  of 
being  able  to  capture  not  just  a  general  non-Maxwellian  velocity  distribution  space,  but  also  radial 
variation  in  mass  flux. 

B.  Hybrid  Potential  Solver 

Generally  two  methods  exist  for  obtaining  plasma  potential  for  plume  simulations.  The  potential  can 
be  obtained  by  solving  the  Poisson’s  equation.  This  approach  requires  cell  spacing  small  enough  to 
capture  gradients  in  the  solution.  Since  the  plume  is  quasineutral,  and  not  magnetized,  a  simpler 
method  is  available  based  on  the  direct  inversion  of  the  Boltzmann  relationship  for  electrons,  cp 
=  cpo  T-  kTe  ln(rLi  /n0).  Although  this  approach  computes  the  correct  potential  in  the  plume 

region,  it  does  not  take  into  account  the  non-  neutral  sheath  region.  In  the  case  of  a  GEO  satellite,  the 
sheath  is  not  negligible  and  can  expand  distances  the  scale  of  the  spacecraft. 


Hence,  we  implemented  a  ”QN  switch”  Poisson  solver  based  on  the  previous  work  of  Santi  and 
Cheng.26  Prior  to  commencing  the  solver  iterations,  the  solver  calculates  the  local  Debye  length  at 
each  node  of  the  simulation  domain.  If  4/3 ttA3  <  V  ,  where  V  is  the  cell  volume,  the  node  is  flagged  as 


quasineutral  and  potential  on  it  is  fixed  to  the  value  obtained  by  the  direct  inversion.  The  Poisson 
solver  then  backfills  the  remaining  region.  Comparison  between  the  two  solutions  is  shown  in  Figure  7. 
The  plots  were  generated  by  injecting  particles  for  a  small  number  of  time  steps  without  using  any 
field  solver.  The  simulation  was  then  restarted  for  zero  number  of  time  steps,  which  resulted  in  the 
potential  update,  but  no  particle  push.  Hence,  the  charge  densities  are  identical.  The  figure  on  the  left 
shows  the  solution  obtained  by  the  inversion  alone.  The  figure  on  the  right  shows  the  solution  from 
the  Poisson  solver  with  QN-switch.  It  should  be  noted  that  a  non-switched  Poisson  solver  was  not  able 
to  converge  for  this  particular  problem.  As  we  can  see,  the  solutions  are  identical  in  the  plume 
region,  as  expected.  The  QN-switch  approach  however  correctly  captures  the  potential  drop  outside 
the  negatively  charged  solar  panel.  The  direct  inversion  method  shown  on  left  effectively  compresses 
the  sheath  to  the  thickness  given  by  a  simulation  cell.  The  ions  are  not  aware  of  the  solar  panel  until 
they  reach  the  cell  adjacent  to  it.  This  difference  has  a  profound  implication  on  the  trajectories  of 
the  charge  exchange  ions,  trajectories  of  which  are  primarily  influenced  by  the  electric  fields  between 
the  plume  and  the  spacecraft  components. 


(a)  Boltzmann  Inversion 


(B)  Poisson  Solver  with  QN  switch 


Figure  7.  Comparison  of  electric  potential  solution  for  an  identical  charge  density  computed  using  the  two 
approaches.  The  solution  is  identical  in  the  plume  region,  but  the  QN-switched  solver  also  resolves  the  sheath 

around  the  negatively  charged  solar  panels. 


C.  Plume  Results 

We  used  this  approach  to  model  the  plume  environment  around  a  generic  spacecraft  operating  a 
cluster  of  two  CHT  thrusters.  The  results  can  be  seen  in  right  section  of  Figure  1.  Prominent 
feature  of  the  spacecraft  is  a  large  solar  array,  which  in  our  model  was  assumed  to  float  negative  in 
respect  to  the  bus.  For  simplicity,  uniform  potential  was  applied  across  the  solar  wing.  The  potential 
on  the  solar  array  was  -20V,  and  the  potential  on  the  bus  was  0V.  The  potential  at  the  thruster  exit 
was  10V.  The  simulation  was  performed  on  a  stretched  40  X  40  X  46  mesh.  The  mesh  extended  from 
8  cm  behind  the  thruster  exit  plane  to  60  cm  in  front  of  the  thruster.  Reference  values  for  the 
potential  solver  corresponded  to  the  potential,  density,  and  temperature  at  the  thruster  exit.  The 
simulation  took  approximately  2  hours  to  complete.  An  actual  detailed  analysis  of  the  plume 
environment  produced  by  this  thruster  is  reserved  for  future  work. 


VII-  Conclusion  and  Future  Work 


In  this  paper  we  presented  a  new  model  for  performing  multiscale  modeling  of  Hall  thrusters.  The 
main  feature  of  our  model  is  its  ability  to  self-consistently  determine  electron  mobility,  plasma 
properties  in  the  thruster  channel,  and  also  the  plume  environment  induced  by  the  thruster  without 
relying  on  supercomputer  resources.  We  demonstrated  the  approach  on  the  2.6cm  Princeton 
Cylindrical  Hall  thruster.  Our  approach  relies  on  a  kinetic  code  that  computes  the  spatial  variation 
of  mobility  by  considering  the  cyclotron  motion  of  electrons,  a  2D  axi-symmetric  code  for  thruster 
discharge,  and  a  3D  plume  code  to  model  the  plasma  environment  and  contamination  effects. 

In  this  paper  we  used  HPHall  to  model  the  thruster  discharge.  As  part  of  our  future  work,  we  plan  to 
perform  the  thruster  modeling  using  our  in-house  designed  2D  code.  Primary  motivation  for  this 
effort  is  to  simplify  the  iterative  processing  needed  to  obtain  a  tmly  self-consistent  solution.  HPHall  is 
written  in  C,  following  non  object-oriented  procedural  implementation.  HPHall  also  contains  physics 
additional  to  what  is  needed  to  capture  the  thruster  plasma  environment.  The  goal  for  our  work  is 
to  create  a  light-  weight  replacement  to  HPHall  that  utilizes  modern  software  engineering  paradigms 
and  easily  ties  in  with  our  kinetic  code. 


Figure  8.  Plasma  density  and  ion  velocity  streamlines  in  the  sheath  with  converging  magnetic  field  lines.  The 
thruster  wall  forms  the  upper  boundary,  and  the  bottom  boundary  extends  into  the  quasineutral  bulk  plasma. 

In  addition,  we  are  investigating  few  additional  components  of  Hall  thruster  discharge.  Modern 
Hall  thrusters,  including  the  CHT,  use  magnetic  fields  that  intersect  the  camber  walls  at  off-normal 
angles.  Configurations  with  highly  oblique  angles  have  been  proposed  to  create  the  so-called 
magnetic  lens,  which  effectively  focuses  ions  and  pushes  them  away  from  walls.  However,  in  such  a 
configuration,  the  radial  component  of  electric  field  normal  to  the  magnetic  field  lines  can  exceed 
the  component  due  to  the  sheath  drop.  In  that  case,  the  net  electric  field  points  away  from  the  walls, 
and  the  sheath  collapses.  We  have  developed  a  code  specifically  to  study  this  near  wall  sheath  region 
in  highly  oblique  magnetic  fields.  A  typical  solution  is  shown  in  8.  The  goal  of  this  work  is  to 
produce  an  algebraic  model  that  can  be  coupled  into  the  thruster  code  to  correctly  compute  the 
sheath  potential  drop. 

References 

1.  Fife,  J.  and  Martinez-Sanchez,  M.,  ’’Two-Dimensional  Hybrid  Particle-In-Cell  (PIC)  Modeling  of 
Hall  Thrusters,”  24th  International  Electric  Propulsion  Conference ,  Moscow,  Russia,  1995,  pp.  1213- 
1224 

2.  Parra,  F.,  Ahedo,  E.,  Fife,  J.,  and  Martinez-Sanchez,  M.,  ”A  two-dimensional  hybrid  model  of  the 
Hall  thmster  discharge, "  Journal  of  applied  physics ,  Vol.  100,  2006,  pp.  023304. 

3.  Hofer,  R.,  Mikellides,  I.,  Katz,  I.,  and  Goebel,  D.,  ’’Wall  sheath  and  electron  mobility  modeling  in 
hybrid-PIC  Hall  thmster  simulations,”  AIAA  Paper ,  Vol.  5267,  2007. 


4.  Koo,  J.  and  Boyd,  I.,  ’’Computational  model  of  a  Hall  thruster,”  Computer  physics  communications , 
Vol.  164,  No.  1-3,  2004,  pp.  442-447. 

5.  Hagelaar,  G.,  Bareilles,  J.,  Garrigues,  L.,  and  Boeuf,  J.,  ’’Two-dimensional  model  of  a  stationary 
plasma  thruster,”  Journal  of  applied  physics ,  Vol.  91,  2002,  pp.  5592. 

6.  Scharfe,  M.,  Thomas,  C.,  Scharfe,  D.,  Gascon,  N.,  Cappelli,  M.,  and  Fernandez,  E.,  ’’Shear-based 
model  for  electron  transport  in  hybrid  Hall  thruster  simulations,”  Plasma  Science,  IEEE  Transactions  on , 
Vol.  36,  No.  5,  2008,  pp.  2058(2068. 

7.  Garrigues,  L.,  Hagelaar,  G.  J.  M.,  Boeuf,  J.  P.,  Raitses,  Y.,  Smirnov,  A.,  and  Fisch,  N.  J.,  ’’Two 
Dimensional  Hybrid  Model  of  a  Miniaturized  Cylindrical  Hall  Thruster,”  30th  International  Electric 
Propulsion  Conference ,  Florence,  Italy,  2007,  IEPC-2007-157. 

8.  Sydorenko,  D.,  Smolyakov,  A.,  Kaganovich,  I.,  and  Raitses,  Y.,  ’’Modification  of  Electron  Velocity 
Distribution  in  Bounded  Plasmas  by  Secondary  Electron  Emission,”  IEEE  Transactions  on  Plasma 
Science ,  Vol.  34,  No.  3,  2006,  pp.  815(824. 

9.  Hirakawa,  M.  and  Arakawa,  Y.,  ’’Particle  simulation  of  plasma  phenomena  in  Hall  thrusters,” 
Proceedings  of  the  24th  International  Electric  Propulsion  Conference ,  Moscow,  IEPC-95-164,  1995. 

10.  Szabo,  J.,  Martinez-Sanchez,  M.,  and  Batishchev,  O.,  ’’Fully  kinetic  Hall  thruster  modeling,”  IEPC 
Paper ,  2001,  pp.  01-341. 

11.  Fox,  J.,  Advances  in  Fully-Kinetic  PIC  Simulations  of  a  Near-Vacuum  Hall  Thruster  and  Other 
Plasma  Systems ,  Ph.D.  thesis,  Massachusetts  Institute  of  Technology,  2007. 

12.  Sullivan,  K.,  PIC  simulation  of  SPT  Hall  thrusters:  high  power  operation  and  wall  effects ,  Ph.D. 
thesis,  Massachusetts  Institute  of  Technology,  2004. 

13.  Raitses,  Y.  and  Fisch,  N.,  ’’Parametric  investigations  of  a  nonconventional  Hall  thruster,”  Physics  of 
Plasmas ,  Vol.  8,  2001,  pp.  2579. 

14.  Brieda,  L.,  Keidar,  M.,  Raitses,  Y.,  and  Fisch,  N.,  ’’Self-Consistent  Calculation  of  Electron  Transport 
in  a  Cylindrical  Hall  Thruster,”  3 1  st  International  Electric  Propulsion  Conference ,  Ann  Arbor,  MI,  2009, 
pp.  1-9. 

15.  Raitses,  Y.,  Smirnov,  A.,  and  Fisch,  N.,  ’’Effects  of  enhanced  cathode  electron  emission  on  Hall 
thruster  operation,”  Physics  of  Plasmas,  Vol.  16,  2009,  pp.  057106. 

16.  Smirnov,  A.,  Raitses,  Y.,  and  Fisch,  N.  J.,  ’’Plasma  measurements  in  a  100W  cylindrical  Hall 
thruster,”  Journal  of  Applied  Physics ,  Vol.  95,  No.  5,  2004,  pp.  2283(2291. 

17.  Brieda,  L.,  Keidar,  M.,  Raitses,  Y.,  and  Fisch,  N.,  ’’Investigation  of  Electron  Transport  in  a 
Cylindrical  Hall  Thruster  using  a  Kinetic  Code,”  45th  AIAA/ASME/SAE/ASEE  Joint  Propulsion 
Conference ,  Denver,  Colorado,  Aug.  2-5,  2009,  2009,  pp.  1-8. 

18.  Wikipedia,  ’’JAVA  Performance,”  http://en.wikipedia.org/wiki/Java_performance 
fCompariso  to  other  languages. 


19.  Particle  In  Cell,  ’’Get  Results  Faster  with  Java  Multithreading,”  http://www. 
particle incell . com/ 2  011/ j  ava-multi threading. 


20.  Lieberman,  M.  and  Lichtenberg,  A.,  Principles  of  plasma  discharges  and  materials  processing , 
Wiley,  2005. 

21.  Nicklass,  A.,  Dolg,  M.,  Stoll,  H.,  and  Preuss,  H.,  ”Ab  initio  energy-adjusted  pseudopotentials  for  the 
noble  gases  Ne  through  Xe:  Calculation  of  atomic  dipole  and  quadrupole  polarizabilities,”  Journal  of 
Chemical  Physics,  1995. 

22.  Szabo  Jr,  J.,  Fully  kinetic  numerical  modeling  of  a  plasma  thruster ,  Ph.D.  thesis,  Massachusetts 
Institute  of  Technology,  2001 

23.  Sydorenko,  D.,  Particle-in-cell  simulations  of  electron  dynamics  in  low  pressure  discharges  with 
magnetic  fields ,  Ph.D.  thesis,  University  of  Saskatchewan,  2006. 

24.  24Brieda,  L.,  Kafafy,  R.,  Pierru,  J.,  and  Wang,  J.,  "Development  of  the  Draco  code  for  modeling 
electric  propulsion  plume  interactions,”  40th  Joint  Propulsion  Conference ,  Fort  Laudardale,  FL,  USA, 
2004 

25.  Nakles,  M.,  Hargus,  W.,  and  VanGilder,  D.,  "Comparison  of  Numerical  and  Experimental  Near-Field 
Ion  Velocity  Distributions  of  the  BHT-200-X3  Hall  Thruster,”  42nd  Joint  Propulsion  Conference , 
Sacramento,  CA,  2006,  AIAA-2006-4479. 

26.  Santi,  M.,  Cheng,  S.,  Celik,  M.,  M.,  M.-S.,  and  J,  P.,  "Further  Development  and  Preliminary  Results 
of  the  AQUILA  Hall  Thruster  Plume  Model,”  39th  Joint  Propulsion  Conference ,  Huntsville,  Alabama, 
2003. 


Section  2,  Plasma-Wall  Interaction  in  Hall  Thrusters  with  Magnetic  Lens 

Configuration 


Introduction 

Hall  thrusters  are  spacecraft  propulsion  devices  that  utilize  applied  magnetic  fields  and  closed 
electron  Hall  drift  to  accelerate  quasi-neutral  plasma.  The  typical  Hall  thruster  consists  of  an  annular  or 
cylindrical  chamber  with  one  end  open  to  the  ambient  environment.  Neutral  propellant  is  injected  through 
the  closed  end.  This  end  also  contains  the  anode.  An  externally  located  cathode  produces  electrons, 
fraction  of  which  enters  the  chamber  and  ionizes  the  propellant.  In  order  to  increase  the  electron  transit 
time,  and  hence  improve  the  ionization  efficiency,  a  magnetic  field  is  applied  over  a  section  of  the 
chamber.  This  magnetic  field  restricts  the  axial  motion  of  electrons,  since  the  electrons  become  trapped  in 
a  closed  azimuthal,  or  Hall,  drift  about  the  thruster  centerline.  The  magnetic  field  thus  also  plays  an 
important  secondary  role.  Since  the  motion  of  electrons  is  restricted  across  the  field  lines,  electrons  will 
tend  to  redistribute  radially  along  the  magnetic  lines  according  to  the  spatial  variation  in  ion  density.  The 
magnetic  field  lines  thus  become  lines  of  constant  potential  and  an  electric  field  develops  in  the  direction 
normal  to  the  magnetic  field.  This  electric  field  accelerates  the  ionized  propellant  out  of  the  device. 

In  the  classical  Hall  thruster,  the  magnetic  field  consists  primarily  of  a  radial  component.  Such  a 
configuration  appears  ideal  at  first  since  it  produces  an  electric  field  with  axial  orientation.  However,  the 
presence  of  walls  modifies  the  near-wall  potential  structure  and  results  in  a  local  component  accelerating 
ions  into  the  walls.  Ion  wall  flux  contributes  to  loss  of  thruster  efficiency  and  to  limited  thruster  lifetime 
due  to  channel  erosion.  In  order  to  mitigate  wall  losses,  some  novel  Hall  thrusters  [27,28]  have  begun 
experimenting  with  magnetic  fields  with  convex  geometry.  Near  the  walls,  this  so-called  magnetic  lens 
induces  an  electric  field  with  a  radial  component  directed  towards  the  channel  centerline  [29].  An 
interesting  aspect  of  the  lens  configuration  is  that  in  the  vicinity  of  the  wall,  the  resulting  magnetic  field 
lines  can  approach  the  wall  with  a  highly  inclined  incidence  angle,  0,  as  measured  from  the  wall  normal. 

Such  a  configuration  generates  an  electric  field  with  a  strong  radial  term  that,  in  the  case  of  a  sufficiently 
large  0,  dominates  the  component  due  to  the  sheath  potential  drop  [30].  This  can  be  seen  from  a  simple 

example.  Consider  a  typical  300V  Hall  thruster  with  a  200V  potential  drop  occurring  across  a  1cm  wide 
acceleration  zone.  The  magnitude  of  the  electric  field  E±  is  then  2xl04  V/m.  Next  consider  the  potential 

drop  due  to  the  wall  sheath.  The  electric  field  along  the  magnetic  field  line  in  the  vicinity  of  the  wall  can 
be  estimated  from  =  Ted\nn / dr~Te/Ar~20eV /0.01cm~2  X  105P/7n  [31].  Here  A r  is  the  sheath 

thickness,  which  is  taken  to  be  10  Debye  lengths.  The  angle  at  which  the  radial  component  of  the  electric 
field  becomes  negative  is  given  by  Er  =  £jj  cos  6  —  E±  sin  0  ~  tan  0  =  E^  /EL,  or  0  ~  85° . 


Ions  are  then  accelerated  away  from  the  wall  and  a  complete  sheath  collapse  is  expected. 
Although  plasma- wall  transition  has  been  the  subject  of  much  past  research,  such  research  typically 
considered  only  the  generalized  radial  case  [32].  In  this  paper  we  investigate  the  sheath  formation  and 
collapse  in  the  presence  of  a  two  dimensional  magnetic  field.  The  analysis  is  performed  using  a  2D 
particle  in  cell  (PIC)  code.  We  use  the  code  to  determine  the  structure  of  the  plasma  sheath  for  several 
magnetic  field  configurations.  We  first  investigate  the  response  of  the  sheath  to  an  inclined  magnetic 
field.  Next  we  extend  the  analysis  to  include  the  influence  of  secondary  electrons  and  a  magnetic  mirror. 


We  also  develop  a  simple  potential  solver  based  on  the  quasineutral  approximation  present  in  standard 
Hall  thruster  codes  to  investigate  the  effect  the  field  solver  has  on  the  sheath  solution.  We  conclude  the 
paper  with  an  analysis  of  sheath  stability,  wall  flux,  and  channel  erosion. 

Computational  Model 

Simulation  Domain 

We  study  the  sheath  formation  using  a  simple  axisymmetric  electrostatic  particle  in  cell  (ES-PIC) 
code.  The  code  is  based  on  the  hybrid  approach  in  which  ions  are  treated  as  particles,  but  electrons  are 
represented  by  a  fluid  model.  The  computational  domain  is  limited  to  a  small  region  near  the  outer  wall, 
as  illustrated  in  Figure  1.  The  inset  illustrates  a  hybrid  annular/cylindrical  Hall  thruster,  in  which  the 
magnetic  field  geometries  of  interest  can  be  found.  The  thruster  schematic  is  based  on  the  Princeton 
Cylindrical  Hall  Thruster  [27].  The  region  being  studied  is  also  highlighted  in  the  inset.  The  small  size  of 
the  domain  allows  us  to  resolve  the  Debye  length  and  thus  directly  compute  the  electric  potential  in  a 
reasonable  amount  of  time  (each  simulation  takes  approximately  20  minutes).  The  domain  captures  the 
acceleration  region  characterized  by  the  presence  of  the  strong  applied  magnetic  field.  In  our  formulation, 
the  anode  and  the  primary  ionization  zone  are  located  to  the  left.  The  upper  boundary  represents  the  wall, 
while  the  bottom  boundary  extends  into  the  quasineutral  bulk  plasma  region.  Ions  are  injected  into  the 
simulation  along  the  left  boundary  and  leave  through  the  open  right  and  bottom  face  or  by  recombining 
with  the  upper  wall. 


Figure  1.  Schematic  of  the  computational  domain.  Ion  particles  are  injected  from  the  left.  The  inset  shows  a 
cylindrical  Hall  thruster  and  highlights  the  region  analyzed  by  our  code. 

To  simplify  the  subsequent  computation,  we  select  a  simulation  mesh  in  which  the  radial 
gridlines  are  aligned  with  the  magnetic  field.  Such  a  formulation  allows  us  to  specify  the  necessary 
reference  values  as  a  function  of  the  axial  grid  coordinate  only.  In  constructing  the  mesh,  we  paid 
attention  to  two  requirements.  First,  the  mesh  had  to  be  capable  of  capturing  the  magnetic  topology  of 
interest:  varying  angle  of  magnetic  field,  and  also  the  magnetic  mirror  effect.  Secondly,  the  mesh  had  to 
be  suitable  from  the  computational  perspective.  Particle  methods  require  scattering  of  particles  to  the  grid 
nodes,  and  conversely  gathering  forces  by  collecting  values  from  the  grid  onto  particle  locations. 
Topologically  structured  meshes  are  preferred  here,  since  physical  coordinates  can  be  mapped  to  the 
computational  space  via  evaluation  of  analytical  functions.  The  mesh  shown  in  Figure  1  satisfies  both  of 
these  requirements.  The  mesh  coordinates  are  given  by 


r  =  r0  +  j  *  Ar 


z  =  i  *  Azj  —  (nr  —  1  —  j)  *  Ar  *  tan(0)  —  0.5  *  (nz  —  1)  *  (A z}  —  A zw) 


(i) 


where  A z}  is  the  local  cell  spacing.  The  cell  spacing  varies  linearly  between  the  top  and  bottom  boundary. 

These  mesh  coordinates  can  be  easily  inverted.  The  j  component  is  obtained  first  from  the  radial  r 
coordinate.  The  i  coordinate  is  then  recovered  from  the  axial  position  z  using  the  second  equation. 

Particle  Injection 

Xenon  ions  are  injected  into  the  simulation  domain  along  the  left  boundary  with  initial  velocity 
uz  =uo  +uth-  Here  u0  is  the  drift  component  and  uth  is  a  random  thermal  velocity  obtained  by  sampling 

the  Maxwellian  distribution  function  at  leV.  The  magnitude  of  the  drift  component  was  set  to  6  km/s, 
corresponding  to  approximately  25eV  of  upstream  acceleration.  Initial  radial  velocity  is  also  obtained  by 
sampling  the  random  thermal  component.  The  number  of  computational  particles  injected  per  time  step  is 
obtained  from  p  =  mAt/w  =  ntuAmAt/w  where  nt  =  5  x  1016  m"3  is  the  injection  ion  density,  and  w  is 

the  macroparticle  weight.  The  weight  was  selected  such  that  each  cell  contained  a  statistically  significant 
number  of  particles.  Particles  are  loaded  with  a  zero  azimuthal  component.  We  assume  that  no  forces  act 
on  ions  in  the  azimuthal  direction  and  hence  the  cylindrical  equations  of  motion  reduce  to  the  Cartesian 
form.  Ion  positions  are  updated  at  each  time  step  according  to  the  Leapfrog  algorithm  by  integrating  the 
Lorentz  force,  F  =  —  eV(p.  The  magnetic  term  is  omitted,  since  in  a  Hall  thruster,  ions  are  not  magnetized. 

Ions  impacting  the  upper  wall  or  leaving  the  computational  domain  are  removed  from  the  simulation. 
Collisions  are  not  included  as  they  generally  play  only  a  minor  role  in  the  sheath. 

Potential  Solver 

Potential  is  computed  by  solving  the  Poisson's  equation,  £0V2(p  =  —e(nt  —ne  —  n5),  with  the 

three  densities  on  the  right  hand  side  corresponding  to  ions,  primary  electrons,  and  secondary  electrons, 
respectively.  Ion  density  nt  is  computed  directly  by  scattering  positions  of  the  kinetic  ions  to  the 

computational  grid.  In  the  frame  of  reference  of  ions,  electrons  respond  instantaneously  to  any 
disturbance.  Electron  motion  is  also  not  impeded  along  a  magnetic  field  line.  The  time-dependent  and 
convective  terms  on  the  left  hand  side  of  the  electron  momentum  equation  then  vanish,  leaving  us  only 
with  the  force  balance  [34], 

d(p  KTe  dne  KTe  dB 

6  dz  ne  dz  B  dz  (2) 


These  terms  correspond  to  the  electric  field,  gas  pressure,  and  magnetic  field  pressures,  respectively.  This 
equation  can  be  easily  integrated  and  subsequently  inverted  to  obtain  an  expression  for  bulk  electron 
density, 


(3) 


This  is  the  well  known  Boltzmann  relationship  modified  by  the  magnetic  field  strength  term.  This 
term  is  seen  to  reduce  the  electron  density  in  regions  of  an  increasing  magnetic  field  -  this  is  the  magnetic 
mirror  effect.  The  standard  Boltzmann  relationship  is  recovered  if  the  magnetic  field  magnitude  remains 
constant  along  the  field  lines.  It  should  be  noted  that  Equation  3  holds  independently  for  each  magnetic 
field  line.  The  three  constants  with  the  0  subscript  correspond  to  the  reference  density,  potential,  and 
magnetic  field  strength.  These  values  are  unique  and  independent  along  each  line.  For  simplicity,  electron 
temperature  is  assumed  to  remain  constant,  kTe  =  10 eV,  and  there  is  no  variation  in  magnetic  field 

strength  in  the  axial  direction,  dB/dz  =  0.  The  reference  density  is  obtained  by  sampling  the  ion  density 

along  the  bottom  edge  of  the  simulation  domain  where  rtl  =ne  =n0.  The  reference  potential  is  assumed 

to  decay  linearly  between  the  left  and  right  boundaries  following  (p0  =  (pL  —  (20,000V/m)(zw  —  zw0). 

The  strength  of  the  magnetic  field  is  computed  from  the  conservation  of  magnetic  flux,  (pm  =  f  B-ds,  or 

BrAz  =  C,  a  constant  value.  Here  A z  is  the  cell  spacing  at  the  corresponding  r  value.  It  should  be  pointed 

out  that,  as  indicated  by  Equation  3,  the  terms  relating  to  the  magnetic  strength  appear  only  as  a  ratio, 
allowing  us  to  select  an  arbitrary  value  for  the  reference  field. 

The  secondary  electron  density  ns  is  obtained  from  V  •  (nu)  =  0.  Density  of  secondary  electrons 

at  the  wall  is  given  by  ns  w  =  sne  w ,  where  s(Te,  0)  is  the  SEE  yield  [31].  In  our  formulation  we  assume 

isotropic  angular  distribution  and  energy  dependence  based  on  the  linear  relationship  given  in  [33], 


(4) 


For  Boron  Nitride,  the  typical  wall  material  in  conventional  Hall  thrusters,  the  coefficients  <j0  and  E1  are 
0.54  and  40,  respectively.  Ep  is  the  energy  of  the  incoming  particle,  measured  in  eV.  Initial  velocity  of  the 
secondary  electrons  is  taken  to  be  usw  =  (2kTw/nmgy/2 .  Energy  conservation  dictates 
u  =  (2qA0/m)1  leading  to 


(5) 


Boundary  Conditions 


Potential  along  the  top  wall  is  fixed  as  <pw  =  <p0  —  A <pw,  where  the  wall  potential  drop  is  given 


by  [31]  as 


A0w  =  reln 


1  -  s(Te,  6) 


(6) 


The  problem  is  closed  by  prescribing  the  normal  electric  field  E±  along  the  left  and  right  boundaries,  and 

zero  tangential  electric  field,  E{1  =  0  on  the  bottom  boundary.  The  electric  field  along  the  left  and  right 

boundaries  is  non-uniform  for  cases  with  a  diverging  magnetic  field  line  topology.  This  can  be  seen  from 
a  simple  observation  of  the  increasing  distance  between  field  lines  as  one  moves  away  from  the  wall.  The 
magnitude  is  obtained  numerically  by  computing  the  normal  distance  d  to  the  next  magnetic  field  line 
(grid  line)  at  each  node.  The  electric  field  is  then  set  from  =  —A (p0/d.  Potential  is  solved  using  the 

finite  volume  method. 


Quasineutral  Solver 


In  addition  to  the  Poisson  solver,  we  have  also  implemented  an  alternative  method  for  obtaining 
the  potential  distribution.  This  approach  was  developed  in  order  to  approximate  the  solution  from  Hall 
thruster  codes  such  as  HPHall  [35,36].  With  the  exception  of  massively  parallel  fully  kinetic  programs, 
codes  developed  to  model  Hall  thrusters  generally  do  not  solve  the  Poisson's  equation.  For  convergence, 
Poisson  solvers  require  mesh  spacing  fine  enough  to  resolve  to  local  Debye  length.  From  the 
computational  perspective,  resolving  the  Debye  length  is  impractical  on  the  spatial  scale  of  a  typical 
thruster.  Instead,  such  codes  rely  on  a  simplified  quasi  one-dimensional  approach  to  solve  electron 
conservation  equations  in  the  direction  normal  to  the  magnetic  field.  The  axial  variation  in  the  reference 
potential,  </>0,  is  obtained  from  the  solution  of  these  equations.  The  radial  variation  in  potential  is  then 


computed  by  assuming  quasineutrality,  ne  =  n t  in  conjunction  with  the  thermalized  potential  model, 
(p  =  (p*  +  kTe/e\n(n/n0  ).  This  approach  is  analogous  to  the  formulation  used  to  derive  the  relationship 


for  bulk  electron  density,  Equation  3.  This  expression  can  be  inverted  to  obtain 


kTe 

<P  =  00 +  — In 


(V) 


where  ne  =nl  =n.  Our  expression  extends  the  quasineutral  formulation  by  taking  into  account  the 
varying  strength  of  the  magnetic  field. 

Implementation 


Simulations  were  performed  on  a  domain  with  50  cells  in  the  axial  and  40  cells  in  the  radial 
direction.  Cell  spacing  was  set  to  ~AD  =  10  _4  m.  The  simulation  time  step  was  adjusted  automatically  by 

the  code  from  its  initial  value  of  3  x  10-9  s  such  that  ion  particles  traveled  no  more  than  0.33  cell  lengths 

per  time  step.  The  thruster  diameter  was  assumed  to  be  6  cm.  The  simulation  domain  was  initially  empty 
of  particles  and  ions  were  injected  into  the  domain  until  steady  state  was  achieved.  Steady  state  was 
characterized  by  zero  net  change  in  particle  counts.  The  simulation  then  continued  for  additional  1000 
time  steps  during  which  results  were  averaged.  The  typical  number  of  computational  particles  at  steady 
state  was  700,000.  Simulation  results,  including  potential,  number  densities,  particle  velocities,  and  wall 
fluxes  were  then  exported.  A  marching  squares  algorithm  was  implemented  in  the  code  to  automatically 
contour  the  resulting  velocity  map  to  obtain  the  sheath  boundary  based.  In  this  work,  we  assumed  that  the 
sheath  boundary  corresponds  to  the  contour  where  the  radial  component  of  velocity  v  =  vB,  the  Bohm 

velocity.  The  code  was  implemented  in  the  Java  programming  language.  Simulations  were  performed  on 
a  Dell  Precision  workstation  with  eight  CPU  cores.  Each  simulation  case  was  launched  as  an  independent 
thread,  and  a  simple  scheduler  was  implemented  to  allow  concurrent  execution  of  the  simulation  cases. 

Results 

Potential  Distribution  at  Uniform  Density 

Often  we  can  obtain  useful  insight  into  the  solution  by  considering  a  simplified  case  that  can  be 
evaluated  in  a  reduced  computation  time.  In  our  case,  we  investigated  the  potential  distribution  that  forms 
in  the  presence  of  a  completely  uniform  plasma.  These  results  are  illustrated  in  Figure  2.  In  all  cases, 
plasma  density  of  5xl016  m'3  was  used.  The  contours  correspond  to  the  lines  of  constant  potential  and  the 
streamtraces  visualize  the  electric  field.  The  classical  Hall  thruster  with  a  solely  radial  magnetic  field  is 
shown  in  Figure  2a).  As  indicated  previously,  this  configuration  results  in  a  primarily  axial  electric  field. 
However,  near  the  wall,  the  wall  potential  drop  modifies  the  electric  field  structure  such  that  the  electric 
field  becomes  oriented  towards  the  wall.  Ions  located  in  this  near  wall  region  are  then  expected  to  be 
accelerated  to  the  wall  and  subsequently  lost  to  the  wall  recombination.  Figure  2b)  illustrates  what 
happens  when  the  magnetic  field  angle  is  increased  to  30°.  Increase  of  the  magnetic  field  incidence  angle 
results  in  a  compression  of  the  region  containing  the  radial  electric  field.  Analogously,  the  critical 
streamtrace  that  delineates  the  near  wall  region  from  the  bulk  acceleration  zone  moves  closer  to  the  wall. 
Ions  located  below  this  line  are  expected  to  be  screened  from  any  wall  effects  and  are  accelerated  in  the 
direction  normal  to  the  magnetic  field.  Ions  above  this  line  are  lost  to  the  wall. 
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Figure  2.  Potential  with  uniform  plasma  density,  a)  6=0°,  Bw/B0=l  b)  6=30°,  Bw/B°=l  c)  6=0°,  Bw/B0=2  d)  6=30°, 

Bw/Bo=2 


Cases  c)  and  d)  show  the  effect  of  the  magnetic  lens.  The  net  angle  of  the  magnetic  field  is  zero  in  case  c), 
however,  the  magnetic  field  strength  decreases  away  from  the  walls.  The  magnetic  field  strength  at  the 
bottom  edge  is  0.5 Bw.  The  solution  is  approximately  antisymmetric.  The  electric  field  is  seen  to  first 
direct  the  ions  away  from  the  wall  in  the  converging  section  of  the  lens.  Ions  are  then  accelerated  back 
towards  the  wall  in  the  diverging  section.  However,  since  the  field  accelerates  the  ions,  the  ion  velocity 
will  be  greater  in  this  half  and  the  diverging  effect  is  expected  to  be  smaller.  The  predicted  result  is  a  net 
acceleration  of  ions  away  from  the  wall.  Case  d)  extends  the  magnetic  mirror  effect  in  c)  by  including  the 
30°  field  inclination  from  b).  This  configuration  approximates  the  magnetic  field  in  thrusters  such  as  the 
Princeton  Cylindrical  Hall  Thruster.  In  the  CHT,  this  configuration  arises  from  the  difference  in  the 
physical  location  of  the  inner  and  outer  magnets.  Figure  2d)  indicates  that  the  effect  of  combined  field 
divergence  and  field  inclination  is  a  reduction  in  the  defocusing  effect  seen  in  c),  while  at  the  same  time 
reducing  ion  wall  flux,  and  retaining  the  net  axial  acceleration  of  ions. 

Sheath  variation  with  magnetic  field  angle 

We  build  on  these  initial  observations  by  performing  series  of  plasma  simulations.  We  first 
investigated  the  effect  of  an  increasing  magnetic  field  angle  in  the  absence  of  a  magnetic  mirror. 
Secondary  electron  emission  was  not  included  in  this  set,  5  =  0  in  Equation  5.  The  magnetic  field  angle  9 

increased  from  0°  (magnetic  field  normal  to  the  wall)  to  85°  (highly  inclined  configuration).  Results  for 
0°,  40  °,  60  0  and  70  0  are  plotted  in  Figure  3.  Ion  densities  are  shown  using  the  contour  plot.  Velocity 
streamlines  are  also  plotted,  as  well  as  the  sheath  boundary.  The  sheath  edge  corresponds  to  the  contour 
where  the  radial  velocity  component  (i.e.  component  normal  to  the  wall)  reaches  the  Bohm  velocity 
%  =  yjkTe/ml.  The  sheath  edge  is  plotted  by  the  solid  red  line.  We  see  that  in  the  case  of  a  zero 

magnetic  field  angle,  the  solution  obtained  by  our  code  is  similar  to  the  well-known  boundary  layer 
problem.  The  sheath  forms  short  distance  from  the  injection  plane  and  continues  to  grow  as  more  ions  are 
accelerated  from  the  bulk  plasma  towards  the  wall.  This  result  is  somewhat  non-physical,  since  in  a  real 
device,  the  sheath  thickness  will  be  finite  at  the  entrance  to  the  acceleration  zone.  Since  our  simulation 
resolves  only  a  small  subset  of  the  Hall  thruster  channel,  we  are  unable  to  capture  the  sheath  that  forms 
upstream  of  our  domain.  To  investigate  the  role  the  initial  sheath  profile  has  on  results,  we  tested  a 
modified  particle  loading  algorithm  in  which  we  allowed  the  injection  density  to  decay  exponentially 
towards  the  wall.  We  found  that  differences  between  the  two  solutions  were  limited  to  several  cells  near 
the  left  boundary.  No  significant  differences  in  plasma  parameters  or  the  sheath  profile  were  observed  at 
downstream  locations.  This  finding  can  be  explained  by  realizing  that  any  ions  injected  into  the  sheath 
will  be  rapidly  lost  to  the  wall.  Hence,  we  ignore  this  initial  region  and  characterize  the  sheath  by  its 
maximum  thickness. 


Figure  3.  Simulation  results  showing  the  sheath  profile  for  three  different  magnetic  field  line  angles,  0°,  40°,  and 
70°,  respectively.  Streamtraces  show  ion  trajectories.  The  red  lines  correspond  to  the  sheath  edge  as  computed  with 

the  Poisson  solver  (solid)  and  the  QN  model  (dashed  line). 

From  our  simulation  we  see  that  the  effect  of  an  increasing  magnetic  field  angle  is  to  reduce  the 
sheath  thickness,  as  stipulated  in  the  previous  section.  At  40°,  the  sheath  thickness,  as  characterized  by  the 
ion  velocity,  asymptotes  to  a  constant  thickness.  The  ion  density  decays  near  the  wall,  however  this  effect 
is  due  to  the  net  acceleration  of  ions  away  from  the  wall  by  the  electric  field,  and  not  by  the  loss  of  ions  to 
the  wall  as  is  the  case  with  the  classical  sheath  model.  At  70°  we  start  seeing  first  evidence  of  a  sheath 
collapse.  The  sheath  has  reduced  in  axial  size  and  extends  only  over  a  small  section  of  the  wall.  In 
addition,  near  the  wall,  the  ion  velocity  streamlines  become  parallel  to  the  wall.  Ions  are  thus  accelerated 
along  the  wall,  instead  of  being  attracted  into  the  wall  as  is  the  case  with  a  smaller  9. 

In  this  study  we  also  looked  at  the  effect  the  potential  solver  plays  on  the  sheath  solution.  The 
sheath  profile  indicated  by  the  dashed  red  line  corresponds  to  the  solution  obtained  using  the  simplified 
quasineutral  approach  similar.  We  can  see  that  the  sheath  thickness  is  artificially  compressed,  and  extends 
to  approximately  three  computational  cells.  For  illustration,  the  computational  mesh  is  displayed  in  the 
top  figure.  The  potential  solution  obtained  in  the  case  of  a  QN  solver  is  completely  driven  by  the  local 
deviation  from  reference  density.  Since  potential  on  each  node  is  evaluated  independently,  the  wall 
potential  is  artificially  screened  at  a  distance  equal  to  one  cell  length.  This  is  in  contrast  with  the  solution 
from  the  Poisson  solver,  which  predicts  similar  screening  to  occur  over  distance  of  several  Debye  lengths. 
The  QN  solution  is  non-physical,  since  the  magnitude  of  the  near-wall  electric  field  is  directly  related  to 
the  cell  spacing.  The  wall  effect  in  the  QN  solver  is  extended  to  the  additional  cells  by  motion  of  ions, 
which  by  their  loss  to  the  wall  reduce  the  local  density.  However,  this  effect  has  only  a  limited  capability 
to  communicate  the  wall  potential  drop  to  the  bulk  plasma.  As  such,  the  QN  solution  results  in  the  wall 
having  a  reduced  ion-attractive  capability  compared  to  the  physically-sound  Poisson  solution.  This  result 
is  demonstrated  in  the  relative  reduction  in  the  sheath  thickness. 


Influence  of  Secondary  Electron  Emission 


We  next  included  secondary  electron  emission  (SEE).  Secondary  electron  emission  is  an 
important  process  in  SPT-type  Hall  thrusters,  in  which  the  acceleration  channel  is  lined  with  an  insulator 
material  for  which  the  secondary  electron  yield  can  approach  unity.  SEE  may  be  an  important  driver  in 
the  so-called  anomalous  electron  transport  across  magnetic  field  lines.  This  effect  is  not  investigated  in 
this  work.  Instead,  we  only  concentrate  on  the  role  of  SEE  on  the  sheath  profile.  Same  set  of  cases 
presented  in  the  previous  paragraph  was  run  with  the  wall  potential  modified  by  the  presence  of  SEE.  We 
include  secondary  electrons  in  our  code  by  computing  the  SEE  emission  coefficient  using  Equation  4. 
From  the  wall  potential  relationship,  Equation  5,  we  can  see  that  the  presence  of  SEE  acts  to  decrease  the 
sheath  potential  and  hence  the  sheath  thickness.  This  prediction  is  confirmed  by  the  result  illustrated  in 
Figure  4.  This  figure  shows  the  potential  contours  for  the  40°  magnetic  field  inclination.  The  dashed  lines 
correspond  to  the  case  with  secondary  electron  emission.  Presence  of  SEE  is  seen  to  reduce  the  sheath 
thickness,  and  thus  is  expected  to  reduce  the  flux  of  ions  to  the  walls. 


Figure  4.  Potential  profile  for  40°.  Dashed  lines  indicate  solution  with  secondary  electron  emission. 

Magnetic  Mirror  Effect 

The  final  set  of  simulations  was  run  to  investigate  the  effect  of  magnetic  mirror.  We  used  mirror 
ratio  Bw  to  B0  of  1.8,  and  compared  case  without  and  with  a  40°  magnetic  field  inclination.  Here  Bw  is  the 
field  strength  at  the  wall,  and  B0  is  the  field  at  the  bottom  edge  of  the  simulation  domain.  We  can  see  from 
Figure  5  that  the  mirror  plays  a  role  similar  to  that  of  the  inclined  field.  As  the  mirror  ratio  is  increased, 
the  divergence  of  the  field  also  increases  and  results  in  acceleration  of  ions  away  from  the  wall.  However, 
unlike  in  the  case  of  a  uniform  magnetic  field,  the  orientation  of  the  electric  field  past  the  centerline  is 
reversed,  and  ions  are  accelerated  towards  the  wall.  This  can  be  seen  in  Figure  5a),  in  which  ions 
trajectories  that  were  initially  parallel,  or  leading  away  from  the  walls,  are  seen  to  turn  towards  the  wall 
past  the  field  centerline.  This  outcome  is  somewhat  different  that  predicted  by  the  case  based  solely  on  the 
magnetic  field  orientation  in  Figure  2c).  We  see  that  the  presence  of  ions,  and  hence  non-uniform  plasma 
density,  results  in  the  sheath  forming  normal  to  the  magnetic  field  lines,  instead  of  parallel  to  the  wall  as 
was  the  case  in  Figure  2.  Figure  5b)  plots  the  solution  obtained  by  including  a  40°  magnetic  field 
inclination.  Ions  are  again  seen  to  be  initially  accelerated  away  from  the  wall.  However,  past  the 
centerline  the  electric  field  becomes  oriented  axially,  acting  to  accelerate  ions  along  the  thruster 
centerline.  This  net  axial  acceleration  reduces  flux  of  ions  to  the  wall  and  reduces  the  sheath  thickness. 


Figure  5.  Ion  density  contours  in  the  presence  of  magnetic  mirror.  Magnetic  mirror  strength  of  2  is  used  in  both 

cases.  Case  (b)  includes  a  40°  magnetic  field  inclination. 


Discussion 

Sheath  Collapse 

In  Figure  6  we  plot  the  variation  in  maximum  sheath  thickness  with  the  incidence  magnetic  field 
angle.  We  consider  three  cases:  inclined  magnetic  field,  inclined  magnetic  field  with  a  magnetic  mirror 
and  SEE,  and  an  inclined  field  solved  using  the  quasineutral  potential  solver.  We  see  that  in  all  cases  the 
sheath  thickness  decreases  as  the  magnetic  field  incidence  angle  is  increased.  Figure  3  shows  that  at  70° 
the  sheath  surrounds  only  a  small  portion  of  the  wall.  From  Figure  5  we  also  see  that  the  sheath  has 
reduced  in  maximum  thickness  by  80%  from  the  value  obtained  at  0°.  The  sheath  thickness  obtained  in 
the  presence  of  magnetic  mirror  closely  follows  the  trend  of  pure  field  inclination.  The  slight  increase  in 
the  thickness  with  the  mirror  may  be  due  to  the  increase  of  the  domain  size  associated  with  the  diverging 
field.  Of  greater  interest  is  the  solution  obtained  using  the  QN  approach.  This  resulted  is  plotted  used  the 
dashed  line.  We  can  see  that  the  QN  method  under-predicts  the  sheath  thickness  by  a  factor  of  3.  In 
addition,  the  QN  approach  also  results  in  a  much  smaller  dependence  of  sheath  thickness  on  the  applied 
magnetic  field  angle.  This  result  indicates  that  the  quasineutral  formulation  used  in  typical  Hall  thruster 
codes  may  not  be  a  good  candidate  for  determine  wall  fluxes,  and  additional  steps  should  be  taken  to 
correct  the  wall  ion  flux[37]. 


Figure  6.  Variation  of  maximum  sheath  thickness  with  the  angle  of  magnetic  field.  Dashed  line  corresponds  to  the 

solution  obtained  assuming  quasineutrality. 


The  sheath  is  seen  to  collapse  at  angles  greater  than  80°,  confirming  the  simple  analysis  presented 
in  the  introduction.  To  better  illustrate  the  dynamics  at  this  highly  inclined  magnetic  field  geometry,  we 
plot  ion  velocity  contours  and  velocity  streamlines  at  the  85°  incidence  angle.  This  resulted  is  shown  in 
Figure  7.  The  contour  plot  corresponds  to  the  radial  velocity  component  normalized  by  the  Bohm  speed. 
We  can  see  that  at  this  high  incidence  angle,  the  radial  component  of  ion  velocity  never  reaches  the  Bohm 
speed.  In  addition,  ions  are  moving  towards  the  wall  only  along  a  small  region  near  the  left  boundary. 
This  result  is  likely  a  direct  byproduct  of  our  loading  scheme  since  it  affects  only  the  ions  injected  into  the 
sheath.  Ions  in  the  bulk  plasma  region  are  instead  accelerated  away  from  the  wall.  Ions  located  just  a 
small  distance  from  the  wall  are  seen  to  follow  trajectory  first  parallel  to  the  wall,  and  subsequently 
turning  away  from  the  wall.  Ions  are  thus  seen  to  be  repelled  by  the  wall,  indicating  a  full  sheath  collapse. 


Figure  7.  Plots  of  normalized  radial  velocity  and  ion  velocity  streamlines  at  0=85°. 

Erosion  and  Lifetime 

Our  numerical  results  confirm  that  the  presence  of  highly  inclined  magnetic  fields  results  in  a 
decreased  flux  of  ions  to  the  wall,  and  eventually  a  sheath  collapse.  This  observation  has  a  profound 
effect  on  both  the  ionization  efficiency  and  the  thruster  lifetime,  since  ion  losses  to  the  walls  are  a  major 
contributor  to  both  of  these  inefficiencies.  Here  we  consider  only  the  wall  erosion.  Material  sputtering 
yield  scales  with  both  the  impact  angle  and  the  energy  of  the  incoming  ions.  Several  models  exist  for 
computing  sputter  yields  for  Boron  Nitride,  the  material  typically  used  in  SPT-type  Hall  thrusters.  In  this 
work  we  utilize  the  logarithmic  fit  suggested  by  Gamier  [38], 

Y0(E)  =  0.0156  \nE  —  0.0638 


This  fit  is  valid  from  the  energy  threshold  of  60  eV  up  to  keV.  In  our  analysis  we  neglect  low  energy 
sputtering.  For  angular  dependence  of  yield,  quadratic  polynomial  fit  is  recommended  by  Yim  [39]. 

Y(E,6 )  =  F0[-4.45  x  1O~704  +  4.91  x  1O'503  -  9.72  x  1O"402  +  3.44  x  1O"30  +  1] 


where  0  is  in  degrees. 

Figure  8  shows  the  computed  wall  flux  and  sputter  yields  for  several  representative  cases.  The 
baseline  is  the  configuration  with  normal  radial  field,  no  SEE,  and  no  magnetic  mirror.  Addition  of 
secondary  electrons  is  seen  to  reduce  the  wall  flux,  although  the  effect  is  small.  Greater  difference  is 
observed  in  the  calculated  sputter  yield.  Since  the  presence  of  SEE  reduces  the  wall  potential  drop, 
impacting  ions  will  possess  lower  energy  and  hence  the  total  erosion  rate  will  decrease  even  if  the  flux 
remained  the  same.  A  much  greater  effect  is  seen  when  the  magnetic  field  angle  is  increased  to  40°.  Flux 
reduces  by  approximately  60%,  leading  to  a  correspondingly  similar  reduction  in  erosion  rate.  Inclusion 


of  the  magnetic  mirror  effect  (at  Bw/B0=  1.8)  reduces  the  flux  even  further,  although  the  effect  is  not  as 
pronounced  as  due  to  the  magnetic  field  inclination. 
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Figure  8.  Comparison  of  wall  flux  and  computed  sputtered  yield  for  5  selected  configurations. 


The  dashed  line  with  markers  corresponds  to  a  40°  incidence  angle  and  the  quasineutral  field 
solver.  This  case  thus  corresponds  to  the  solution  plotted  with  the  dotted  line,  SEE  and  40°.  We  can  see 
that  although  the  QN  solver  produces  trend  qualitatively  comparable  to  the  reference  case,  the  computed 
fluxes  and  yields  do  not  agree  quantitatively.  Average  QN  flux  is  3.5xl018  #/m2/s  compared  to  7.2*  1018 
#/m2/s  obtained  using  the  Poisson  solver.  This  difference  corresponds  to  a  50%  numerical  reduction  in 
wall  flux.  Comparably,  the  average  sputter  yield  decreases  from  1.5><1017  #/m2/s  to  3.1xl016  #/m2/s  with 
the  QN  solver,  corresponding  to  an  80%  decrease.  The  QN  method  for  obtaining  potential  is  seen  to  both 
decrease  the  wall  flux,  and  also  decrease  the  energy  of  the  impacting  ions,  underpredicting  the  erosion 
rates  by  almost  one  order  of  magnitude. 

Sheath  Stability 

The  analysis  presented  in  the  previous  paragraphs  was  performed  using  the  prescribed  normal 
component  of  electric  field  EL  =  20  kV/m.  The  electric  field  magnitude  is  a  design  parameter  of  the 

thruster,  as  it  is  arises  from  the  particular  design  of  the  magnetic  circuit  and  also  the  potential  drop 
between  the  cathode  and  the  anode.  To  investigate  the  effect  the  field  strength  has  on  the  sheath  profile, 
we  ran  the  code  for  several  values  of  E±  with  0=60°,  no  SEE,  and  no  magnetic  mirror.  These  results  are 

plotted  in  Figure  9.  The  solid  line  at  20kV/m  corresponds  the  case  shown  previously.  Reducing  the 
applied  potential  drop  leads  to  a  thicker  sheath,  as  expected.  We  can  see  that  for  this  particular  case,  full 
sheath  collapse  will  occur  approximately  at  EL  =  50  kV/m.  It  should  be  pointed  out  that  this  model 


predicts  that  the  inclined  magnetic  field  leads  to  sheath  formation  only  if  axial  electric  field  is  small. 
According  to  this  model  it  is  predicted  that  the  potential  drop  inside  the  Hall  thruster  channel  decreases 
leading  to  shifting  the  potential  drop  outside.  Such  effect  has  significant  implications  on  the  plume 
formation  and  the  thruster  contamination  aspects.  It  is  interesting  to  point  out  that  such  trend  is  also 
observed  experimentally [40,4 1 ,42] . 


Figure  9.  Sheath  profile  for  6=60°  as  a  function  of  normal  electric  field 


Conclusion 

In  this  paper  we  investigated  the  topology  of  the  plasma  sheath  that  forms  along  the  wall  in  the 
acceleration  zone  of  a  Hall  thruster  in  the  presence  of  a  two  dimensional  magnetic  field.  Our  analysis 
concentrated  on  the  effect  an  increasing  angle  of  magnetic  field  incidence  plays  on  the  maximum  sheath 
thickness.  It  was  shown  that  at  highly  inclined  angles,  ions  are  repelled  by  the  wall  and  the  sheath 
collapses.  In  addition,  we  studied  the  effect  of  secondary  electron  emission  and  magnetic  mirror.  Both  of 
these  processes  reduced  wall  fluxes,  however,  the  effect  was  negligible  compared  to  the  effect  obtained 
purely  by  the  increased  magnetic  field  incidence  angle.  We  used  the  computed  wall  fluxes  along  with  an 
erosion  fit  to  predict  the  effect  the  magnetic  field  topology  has  on  thruster  lifetime.  The  erosion  rate  was 
seen  to  decrease  by  almost  one  order  of  magnitude  with  increase  of  the  incident  field  to  40°.  In  addition, 
we  compared  our  results  to  those  obtained  by  typical  Hall  thruster  codes  that  utilize  simplified 
quasineutral  approach  to  obtain  the  radial  electric  field.  We  found  that  the  QN  approach  is  not  suitable  for 
determining  wall  fluxes  and  erosion  rates,  as  the  computed  results  underpredicted  the  results  obtained  by 
the  more  physically  accurate  Poisson  solver  by  almost  one  order  of  magnitude. 
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